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ABSTRACT 

Detecting and characterizing the anisotropy pattern of the arrival directions of the highest energy cosmic rays are crucial steps to- 
wards the identification of their sources. We discuss a possible distortion of the cosmic ray flux induced by the anisotropic and 
inhomogeneous distribution of extragalactic magnetic fields in cases where sources of ultrahigh energy cosmic rays are rare transient 
phenomena, such as gamma-ray bursts and/or newly born magnetars. This distortion does not involve an angular deflection but the 
modulation of the flux related to the probability of seeing the source on an experiment lifetime. To quantify this distortion, we con- 
struct sky maps of the arrival directions of these highest energy cosmic rays for various magnetic field configurations and appeal to 
statistical tests proposed in the literature. We conclude that this distortion cannot affect present experiments but should be considered 
when performing anisotropy studies with future large-scale experiments that record as many as hundreds of events above 6 x 10" eV. 
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1. Introduction 

The origin of ultrahigh energy cosmic rays (UHECR) - par- 
ticles with energy E > 10''' eV - remains one of the 
longest running outstanding problems of astrophysics. The 
extremely low flux that prevails at the highest energies, ~ 
1 /km" /century, makes it difficult to accumulate more than a 
handful of events per year at ~ 10^" eV {- 100 EeV). The 
construction of large-scale detectors has nevertheless brought 
a wealth of significant results in t he past few years. In par- 
ticular, the Fly's Eye experiment (Abbasietal. 2008, '2010a|) 
and the Pierre Auger Observato ry ( Abraham et al. 2008a; 
iPierre AUGER Collaboration et aLiuOiol) have detected a high 
energy cut-off at ~ 6 x 10'^ eV, in excellent agreem ent with 
the predicted Greisen-Zatsepin-K uzmin (GZK) cut-off dGreisenl 
119661: IZatsepin & Kuzminlll966h . which results from the pion 
production interaction of protons above that energy with cosmic 
microwave background photons. Strictly speaking, one cannot 
exclude at this stage that this cut-off represents the maximal en- 
ergy during acceleration, yet if one accepts the GZK interpreta- 
tion, this result would confirm that the sources of ultrahigh en- 
ergy cosmic rays are extragalactic and located within the GZK 
sphere of radius ~ lOOMpc. Similar arguments and conclusions 
can be drawn for ultrahigh energy nuclei, although the energy 
loss differs, being caused by the photodissociation of the nuclei. 

The Pierre Auger Observatory has also announced the de- 
tection of anisotropy at the 99% c.l. in the arrival direction 
map o f UHECR with £ > 6 x lO^'^eV ([Abraham etal.ll2007l 
l2008bl) . This anisotropy is expected if the (unknown) sources 
of ultrahigh energy cosmic rays reside in large-scale structure 
- which itself appears anisotropic up to depths of several hun- 
dreds of Mpc - provided that the magnetic deflection accumu- 



lated by the UHECR is not too large. A detailed analysis of 
the anisotropy pattern indicates that it is consistent with a distri- 
bution of UHECR sources across the local large-scale structure 
djCashti & Waxmanll2008l:lKoers & Tinvakovl2009t[Aublin et al.l 
20Q1). 

If confirmed by future data - this anisotropy is at present 
rejected by the HiRes experiment in the northern hemisphere 
( Abbasi et al. 2010bJ) - this result would be of prime importance, 
notably because it would open the way to charged particle as- 
tronomy. In particular, one could envisage probing the distribu- 
tion of the population of ultrahigh energy cosmic ray sources 
with future large-scale cosmic ray detectors, a crucial step to- 
wards the identification of the sources of these particles. This, of 
course, would require a large effort on the experimental side, to 
build detectors with an aperture significantly larger than already 
achieved. On the theoretical side, it also requires the understand- 
ing of the various biases that can distort our reconstruction of 
source properties from sky maps of arrival directions. 

Magnetic fields on large scales. Galactic and intergalactic, 
are prime suspects in this regard, and many studies have been 
devoted to characterizi ng their influence on the ar rival directions 
(see the discussion in iKotera & Lemoind l2008l and references 
therein). As our understanding of cosmic magnetic fields has 
improved over the years, the models have evolved. In particular, 
it is now understood that the typical amount of deflection must 
vary according to the direction in the sky, not only because the 
Galactic mag netic field exerts a different influence in differen t 
directions (eg lAlvarez-Muniz et al ] l2001llTakami & Satdl2008h . 
but also because the projected extragalactic magnetic field ap- 
pears to be anisotropic, as does the large-scale structure on the 
GZK sphere length scale. 
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In the present paper, we discuss a possible distortion of the 
maps of arrival directions that is associated with the inhomoge- 
neous distribution of large-scale magnetic fields which has so 
far not been discussed. If sources of ultrahigh energy cosmic 
rays are rare and powe r ful events, such as gamma- ray bursts 
dMilgrom & Usovlll995l: I Vietril 1 1 9951; I WaxmanI 1 1 9951) or newly 
born extragalactic magnetars ( Aronsll2003h . the probability of 
detection of cosmic rays from a given source is extremely small 
in an experiment lifetime unless the arrival times of these cosmic 
rays have been substantially dispersed by intervening intergalac- 
tic magnetic fields. As we discuss in Section 2.1, this implies 
that the flux in those directions of the sky with small integrated 
column densities of magnetized gas should be low relative to the 
mean background flux. To assess the magnitude of this effect, 
we perform numerical Monte Carlo simulations of the event de- 
tection, assuming that the source distribution of ultrahigh cos- 
mic rays f ollows the large-scale structure (traced using the PSCz 
catalogue. ISaunders et al.ll2000h and accounting for the inhomo- 
geneous distribution of intervening magnetic fields. 

Our study indicates that the above-mentioned effect could 
only marginally affect the results of current experiments such 
as the Pierre Auger Observatory, but that it should be taken 
into account by future experiments that might record as many 
as 10^ - lO-' events above 60EeV. The magnitude of this effect 
is directly controlled by the average depth to scattering against 
magnetized structures in the intergalactic medium. The layout of 
the present paper is as follows. In Section 2, we discuss the gen- 
eral conditions relating to the source characteristics, the degree 
of magnetization of the Universe for which the above effect ap- 
plies and we specify how the bias is modeled. In Section 3, we 
construct maps of arrival directions of ultrahigh energy cosmic 
rays for various magnetic configurations and we present quanti- 
tative estimates of the magnitude of the effect for simulated sets 
of events. We draw our conclusions in Section 4. 



2. Model 



2.1. Magnetic field topology 



As discussed in detail in iKotera & Lemoin3 (l2008f) . the distri- 
bution of magnetic fields in the Universe, as experienced by a 
ultrahigh energy cosmic ray, is quite likely to be highly inho- 
mogeneous, being comprised of magnetized structures such as 
filaments, halos of starburst galaxies, lobes of radio-galaxies, 
and halos of clusters, separated by giant unmagnetized voids. As 
long as the magnetic field B < 10 '^G in some region, this re- 
gion can be considered as unmagnetized from the point of view 
of the cosmic rays that traverse it. We note that even extreme 
models that predict magnetogenesis at high redshift with a strong 
magnetic field lead to the above picture, of a Universe in which 
filaments are much more strongly magnetized than the voids, as a 
result of the amplification of the magnetic field during structure 
formation and the dilution of the magnetic field in expanding 
voids. 

In this situation, one may describe the propagation of the 
ultrahigh energy cosmic rays as a sequence of stochastic inter- 
actions with magnetized structures, which play the role of scat- 
tering agents. Notwithstanding deflection by the Galactic mag- 
netic field, the overall effect is characterized in particular by 
the typical optical depth for magnetic deflection, r. The opti- 
cal depth is dominated by the structures with the largest ncr, 
where n represents the space density and cr the cross-section 
of the magnetized halo. Sca nning through the possible scatter- 
ing agents in the Universe, iKotera & Lemoina (l2008h pointed 



out that the major protagonists in this stochastic UHECR prop- 
agation model are: (1) halos of old radio-galaxies with num- 
ber density, typical radius and magnetic field in the range of 
«rg ~ 3 X 10-^ - 3 X 10-2Mpc-^ ~ 1-3 Mpc, and 
fire ~ 1 - 10 X lO ^^G respec tively (iMedina-Tanco & EnBlinI 
120011: iFurlanetto & Loebll200TI) : (2) magnetized galactic winds 
with «gw = 2 - 5 X 10"^ Mpc"^, rgw = 0.5-1 Mpc, and 
Brg ~ 1 - 10 X lO^^G dBertone et al.ll2006l) : and/or (3) the fil- 
aments/walls themselves if they have been magnetized by some 
process or the filling factor of the above radio relics and super- 
winds in the filaments reaches unity. This can be summarized by 
saying that magnetized scattering centers are distributed in fil- 
aments/walls of large-scale structure - with typical interaction 
lengths of df ~ 30 Mpc - and that these magnetized structures 
account for a covering fraction rj < \ of the filament/wall pro- 
jected area. The typical optical depth for magnetic scattering on 
a GZK distance scale then reads 



<r) - 3/7 



/ 



100 Mpc /\ 30 Mpc 



df 



(1) 



It is important to recall that one might well find rj to be sig- 
nificantly smaller than unity if the filling factor of magnetized 
regions is small, as suggested for instance by the latest numer- 
ical simulations of enrichment of the intergalactic medium re- 
ported in Donnert et al. (2009!). Nevertheless, one can bound the 
average optical depth below by the contribution of clusters of 
galaxies to the scattering depth: with a density «c ~ 10"^ Mpc^^*, 
and a radius out to which magnetization can be considered suf- 
ficient for UHECR deflection, rc ~ a few Mpc, one finds that 
<T)e ~ 0.08(r/5Mpc)2(//100Mpc). We use a lower bound of 
0. 1 for (t) to a depth of 100 Mpc in what follows. 

Since the source distance of protons of energy > 60 EeV can- 
not exceed 200 Mpc - comparable or smaller horizons are found 
for heavy nuclei - one expects the optical depth r to that distance 
to depend on direction angle, just as the projected large-scale 
structure. 



2.2. UHECR bursting sources and Inhomogeneous cosmic 
magnetic fields 

Henceforth, we assume that the sources of ultrahigh energy cos- 
mic rays are of the bursting type, with a rate per unit volume 
of «s. If the sources were steady emitters, then the effect de- 
scribed in this paper would not take place. The time delay im- 
parted by the magnetic field cannot modify the flux predic- 
tions of steady emitters; the only effect of intervening mag- 
netic fields is in this case limited to deflection. Nevertheless, we 
note that there are phenomenological arguments against the ex- 
istence of steady sources o f ultrahigh energy cosmic rays (see 
iLemoine & Waxrrianll2009l for instance). 

We consider a volume element 6V located at a distance 
r from us, with an optical depth t against magnetic scatter- 
ing. Cosmic rays emitted by sources in 6V have a probability 
Pt = I - exp(-T) of suffering at least one interaction with a 
magnetized system, hence of suffering a dispersion in their 
arrival times, which can be calculated as follows. One inter- 
action with a structure of magnetic field strength B, coherence 
length Ab, and typical transverse size r, implies a time delay of 
dt ^ 10^yr(r/2Mpc)^(fl /10nG)^(^fi/100 kpc)(£/ 100 EeV)-^ 
(IKotera & Lemoinel2008h : the total time delay then scales as the 
number of interactions A^int (which is Poisson distributed with 
a mean of t). As pion production occurs in a stochastic way, 
the dispersion in the arrival times at energies above 60 EeV is 
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comparable to the total time delay (Waxman & Miralda-Escude 
1996, Lemoine et al. 1997, Kotera & Lemoine 2008), hence 
2f ~ NintSt. The average dispersion for those particles that ex- 
perience at least one interaction can be straightforwardly eval- 
uated as S, - T [1 - exp(-T)]"' 6t. Consequently, a detector of 
area Aexp and exposure time Texp (with Texp ^ cr, expected) will 
record a mean number of A^ev/s,i:, - A^UHECR^exprexp/(47rr-Sf) 
from each source, assuming each produces A^uhecr cosmic rays 
above some threshold energy Ethr- The mean number of sources 
contributing is NsXi - ns^V'Lt, so that the mean total number 
of events received is A^obs.tr, = A^UHECRAexp7'expns(^W(47rr-). As 
expected, the dispersion in arrival times does not influence the 
flux. As is well-known, the cosmic -ray energy output Euhecr 
per source and the source density rate are related to each other, 
in this case, through the normalization of the experimentally de- 
termined cosmic ray flux: ^uhecr'^s - 0.5 x 10"'*'* erg/Mpc^/yr 
dKatz et al.ll2009l for a recent estimate). 

In principle, cosmic rays emitted by these sources also have 
a probability I - pj = exp(-T) of traveling to Earth without suf- 
fering an interaction with a magnetized system in the intergalac- 
tic medium. However, from the point of view of an experiment, 
with lifetime T, the detection of these particles is extremely un- 
likely if the sources are rare and bursting, such as gamma-ray 
bursts. For instance, the occurrence rate of sources within a fi- 
nite solid angle AQ, out to some distance r, is given by: v(AQ) - 
3 X 10"'*yr"i AQr^(),)«s,_9, where = /iJlO"'' Mpc^^yr"', 

'"100 = r/lOOMpc, and AQ is expressed in steradians. For refer- 
ence, the most recent estimate of the local l ong gamma-ray burst 
rate i s 1.3 + 0.6 - 0.7 x 10 ''Mpc^^yr ' dWanderman & PiranI 
120091) . The probability of detecting cosmic rays in that area of 
the sky is then rv(AQ) « 1. 

If such cosmic rays were detected, they would produce a 
giant flare in the detector, i.e., a localized anomalously large 
number of cosmic ray events detected during a short period 
in time, if is as low as that considered above. In the ab- 
sence of intergalactic deflection, the total dispersion in ar- 
rival times indeed reduces to the contribution of the turbu- 
lent component of the Galactic magneti c field; with a scale 
height of ~ 2 kpc dHan et alj l2006h , strength 6Bmw ~ 
3 fiG and coherence length Ab ~ lOOpc, one finds a time de- 
lay 6tMw - 0.2yr(rB/2kpc)2(5BMw/3yuG)-(/lB/100pc)£'^(J and 
o"r,MW ^ i^^MW, assuming that these cosmic rays are protons (see 
further below for heavier nuclei). Cosmic rays can be expected 
to have a similar dispersion in arrival times upon exiting the host 
galaxy of their source. Assuming that the total time dispersion 
for these particles is 6tMw ^ T, the mean number of events 
per source is A^ev/s.o = A^UHECR^exp/(47r'"~), much larger than ex- 
pected according to the previous limit by a factor of S,/rexp. As 
to the mean number of sources contributing to the flux, one finds 
A^s.o - 'h^VTsxp- To provide quantitative estimates, we write 
^UHECR - 10^^ £53 erg (£53 ~ 1 provides the correct normal- 
ization of the flux for «s ~ 10"^ Mpc^-'yr"') and estimate that 
A^ev/s,o ~ 3 X 10'*£'53r^Qg (for Aexp = 3000 km- corresponding to 
the Pierre Auger Observatory). No such flare has been observed 
by the Pierre Auger Observatory, for which the typical number 
of events per source ab ove 60 EeV is on the order of unity or less 
dAbraham et al.ll2008ah . 

When such flares are not detected, and one assumes that the 
sources are rare and bursting, then the arrival direction flux must 
be modulated by the probability of suffering at least one interac- 
tion with a magnetized system, which is pj - 1 -exp(-r) as dis- 
cussed above. We note that this modulation might be even more 
extreme if the particle had to encounter more than one scattering 



event in order to disperse its arrival time sufficiently for it to be 
detectable; we adopt the above (slightly conservative) modula- 
tion in the following. 

The above modulation formally corresponds to the condi- 
tional probability of detection given that no flare has been de- 
tected (under our model assumption that the source is bursting, 
of rare occurrence). It can also be seen as the consequence of a 
hierarchy of timescales. At a continuous rate, the detector regis- 
ters particles that have experienced at least one interaction with 
magnetized systems in the intergalactic medium, provided that 
the occurrence rate v of bursting sources in the GZK sphere 
and the average dispersion across that distance scale satisfy 
X,v » 1 . However, the arrival direction map of these particles 
is modulated by pr- To detect flares of particles that have trav- 
eled without interacting with magnetized structures, one would 
need to integrate over a timescale V"', but V"' » Texp for rare 
bursting sources. 

Obviously, if the energy output of each source is small 
enough (see below), i.e. if the source occurence rate is suffi- 
ciently high, the flares detected without intergalactic deflection 
could produce fewer than one event in the detector, in which 
case one could not tell whether intergalactic deflection has taken 
place, and the modulation 1 - exp(-r) should not be applied. 
The latter statement depends of course directly on the expo- 
sure of the instrument. For the Pierre Auger Observat ory, for 
instance, if no significant clustering is found (Abrah am et al.l 
l2008ah then no such flare (i.e. without intergalactic time dis- 
persion) has been detected from bursting sources with energy 
output /iuHECR ^ 10'*^ erg. This implies that for sources with an 
occurrence rate «: 10"'' Mpc^^yr"', the above modulation 
1 - exp(-T) should be taken into account when comparing the 
arrival directions of ultrahigh energy cosmic rays with possible 
source distributions. For future instruments with an aperture ten 
times that of the Pierre Auger Observatory, this upper bound on 
the density rate should be increased tenfold. 

To summarize, the distortion of arrival directions maps that 
we examine in the following Section takes place for bursting 
sources that rarely occur in the GZK sphere, as measured com- 
paratively to an experiment lifetime. Gamma-ray bursts sources 
of ultrahigh energy cosmic rays and magnetars fall into this 
category. Regarding the latter, the typical energy output is ~ 
10^' ergs at 10^*^ eV for a typical ra te of 10"Vyr per galaxy, 
i.e. ~ 10"^Mpc"^yr"' dAronsll2003h . The acceleration of ul- 
trahigh energy cosmic rays in blazar flares has also been dis- 
cussed dPermer et al.l 120091) ; however, in this scenario the par- 
ticles experience a substantial time delay ~ 10^ yr upon exiting 
the lobes of the associated radio-galaxy, hence the source ap- 
pears as a steady emitter of ultrahigh cosmic ray sources for all 
phenomenological purposes. The effect that we discuss therefore 
does not apply to these models. 

One cannot exclude that ongoing or future detectors could 
detect a giant flare associated with the observation of events that 
have not suffered dispersion in their arrival times, all the more 
so at source occurrence rates of «s ~ 10"^ Mpc"^yr"' and small 
depth to magnetic scattering of r < 1. The typical rise and decay 
timescale then corresponds to < dtMw, a fraction of a year In 
this case, one would be able to asses directly the cosmic ray 
energy output of the source, along with the possible detection of 
a counterpart, which would provide unvaluable information. 

Finally, we should mention that we have assumed that ul- 
trahigh energy cosmic rays are protons. There is currently 
no clear determination of the chemical composition of ultra- 
high energy cosmic rays. While the HiRes experiment has 
reported a composition with a proton fraction that increases 
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beyond the ankle (lAbbasi et al l l2010ah . the measurements of 
the Pierre Auger Observatory point instead toward a com- 
posi tion that becomes increasingly heavier beyond the an- 
kle (iPierre AUGER Collaboration et al.ll2010l) . By themselves, 
these results of the Pierre Auger Observatory are somewhat 
puzzling because they may be in apparent contradiction with 
the anisotropy results (see e.g. iLemoine & Waxmanll2009l) and 
no hadronic model of shower reconstruction has so far been 
able to reproduce a ll of the composition data of the Pierr e 
Auger Observatory dPierre AUGER Collaboration et aP 1201 Oh . 
Nevertheless, if the composition were predominantly heavy with 
charge Z » 1, then the above discussion would remain valid, 
but the sky maps discussed below would not be realistic owing 
to the large deflection imparted by the Galactic magnetic field. 
In practice, our imperfect knowledge of the Galactic magnetic 
field would limit the impact of the distortion discussed below. 

2.3. Flux calculation 

To model the flux and simulate the events, we adapt the method 
introduced bv iWax man et al. (1997): in a given volume ele- 
ment SV, the mean number of cosmic ray sources is written 
NsXi - nJL,6V, where 2, - r^f/ [1 - exp(-r)] on average, 
following section 12.21 in our simulations, we neglect the ran- 
dom nature of S,, which can be regarded as a higher order ef- 
fect as its magnitude does not modify the average flux, only 
the random fluctuations around this average. The mean number 
of cosmic ray events emitted by each source is then, following 
the earlier discussion, A^ev/s,i;, = A^UHECR^exp7'exp/(47rr^2;,) and 
the mean contribution to the flux from this region is therefore 
A^ev/s,s,A^s.i;,- To account for multiple events per source, one con- 
siders the probability that S sources produce i events; these S , 
variables are indeed independent of each other and Poisson dis- 
tributed with mean Nsx,Pi^ where P, - N'^^^^^ exp(-A^ev/s,s,) /'! 
is the Poisson probab ility that a given source produces / events 
(IWaxman et al.lll997h . 

One can use the above direct averages multiplied by 1 - 
exp(-T) to compute an average flux distribution, as we do in 
Section 2.3.1. Direction-dependent effects are incorporated as 
follows. We use the following law to express the optical depth at 
distance r in direction n 



T(r,n) = (r) f d/^^f^, (2) 
Jo <Pg) 

which amounts to scaling the density of magnetized structures 
according to the galaxy density pg(Z, n) (at distance I, in the 
direction n). The average optical depth against scattering with 
magnetized structures (t), is calculated by averaging over all 
sky, and bears the typical values calculated in Eq. ([TJ. In what 
follows, we use the definition (t)ioo to express this quantity 
calculated to a depth of 100 Mpc. The galaxy density is traced 
through the PSCz survey of galaxies (Saunders et al. 2000). We 
also assume that the density of sources follow the galaxy density 
as traced by the PSCz survey. These choices are motivated by the 
PSCz being a good tracer of star-forming galaxies, which should 
host bursting sources of ultrahigh energy cosmic rays such as 
gamma-ray bursts and magnetars, and contribute to the magnetic 
enrichment of the intergalactic medium. We then integrate the 
contributions of each volume element along r in each direction 
of the sky, modulated by the appropriate 1 - exp(-T) probability 
of intersection with magnetized structures, to compute the map 
of arrival directions. 

We note that, for a given average source density «s, the final 
average number of events detected depends on the average opti- 



cal depth as a consequence of the flux modulation by 1 -exp(-T) 
in each direction. To produce comparable sets of events, we thus 
rescale the density as a function of t in order to keep the av- 
erage number of events fixed as t varies. It is straightforward to 
see that the all-sky average flux out to a distance d scales with 
average magnetic depth t as [1 - 1/{t) + exp(-{T))/{T)] (with 
(r) calculated up to depth d) if one assumes a homogeneous and 
isotropic distribution of sources and magnetized structures; if 
(t) « 1, the flux is thus suppressed by (r)/2, but the suppression 
factor converges to unity as 1 - 1 /(r) for (t) > 1 . Therefore, for 
an average optical depth of order unity on the horizon distance 
scale for particles above the threshold energy E±i-, the source 
density must be renormalized by a factor of order unity. This 
renormalization lies well within the margins of the uncertainties 
in the absolute rate of possible sources. 

To simulate a finite sets of events, we resort to Monte Carlo 
simulations, in which we rely on the above probability laws for 
5,: in each volume element of the sky, we draw the random vari- 
able Si, which characterizes the number of sources contributing 
i events. The total number of events detected from that volume 
element is then iS i. In order to account for the modulation 
effect by pr, we assume that the mean number of sources con- 
tributing to the flux in that volume element is PrNsx,- 

2.3.1. Mean flux 

The cosmic ray (CR) differential flux produced by sources at 
distances smaller than rmax in direction n is given by 



N_ 
4^ 



Jo 



dl «s b(l, n) , 



(3) 



where is the injection rate of particles by each source, b de- 
notes the bias of the source distribution, which we relate - as ex- 
plained above - to the galaxy density, i.e. b{l,n) = pg(/, n)/(p), 
(p) is the mean density of the PSCz catalog, «s - '^I'ts is the av- 
erage apparent source density, and 2, is the typical dispersion in 
arrival times (which we keep its average value fixed to, as dis- 
cussed above). We consider two values of the apparent source 
density: - 10^^'Mpc"-' and «s - lO^^Mpc"^. The former cor- 
responds to the number density of bright galaxies and the second 
is a lower bound inferred from experimental data, as follows. 
First of all, the distance traveled by the highest ener gy event 
(E > lO^^eV) detected by the Fly's Eye experiment (iBird et al.l 
1 19941) cannot exceed 50 Mpc; its detection thus indicates that 
at least one CR source exists in the field of view of Fly's Eye 
out to 50 Mpc, implying a source density «s Z lO^^Mpc"-*. 
Furthermore, the non-dete ction of significant clus tering by the 
Pierre Auger Observatory d Abraham et al.ll2p08bl) also implies 
that W s > lO^^Mpc"^ (see the discussion in Kashti & WaxmanI 
I2OO8I) . 

These two different apparent source densities lead to dif- 
ferent predictions of the detection of multiple ts of events. For 
/Is = lO^^Mpc""* for instance, there are fewer sources in the 
GZK sphere, hence for a projected sample of A^ev - 100 par- 
ticles of energy above 80 EeV from sources located within a 
distance below rmax = lOOMpc, the number of apparent sources 
is ~ 40 leading to the appearance of multiplets, which can be 
used to constrain the apparent source density. If ris - lO^^Mpc"^ 
however, the number of sources is ~ 4 x 10^, hence the average 
number of event per source is smaller than unity. For an energy 
threshold of 60 EeV, the horizon increases to 200 Mpc, hence 
for a sample of 1000 particles as considered further below, the 
conclusions are similar. 
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We expect to detect more data from future experiments, such 
as Auger North and JEM-EUSO. The development of Auger 
North is due to start in 201 1 in Colorado (USA); this experiment 
should reach an area of more than 20000 km^, which represents 
an increase of about nearly an order of magnitude when com- 
pared to Auger South jBliimer & the Pierre Auger Collaboration! 
120101) . The Extreme Universe Space Observatory onboard the 
Japanese Experiment Module (JEM-EUSO) is to be developed 
on th e International Space Station in 2014 ( lEbisuzaki et alj 
I2OIOI) . With a field of view of 60° from a 430 km orbiting al- 
titude, the integrated exposure time should exceed 10^ km- sryr. 
One should thus expect to detect more than 10^ particles above 
7 X 10'^ eV during its five year operation. For this experiment, 
the non-observation of multiplets of events of energy above 80 
EeV would indicate that «s > lO^-'Mpc"^. 

In our simulations, the flux is normalized to the expected 
total flux above a threshold energy E, i.e. A^tot(> E) - 
J '4exp7^exp/(f^)^(n)dQ, where as before, Aexp denotes the ef- 
fective detector area, Texp the observation time, and f(Q.) the 
aperture of our fiducial experiment. 



3. Simulations 

3.1. Sky maps 

Using Eq. (O, we construct sky maps of the average UHECR 
flux for two energy thresholds: 60 EeV and 80 EeV, correspond- 
ing to integration depths of 200 Mpc and 100 Mpc, respectively. 
These sky maps, that do not include any experiment exposure 
factor for completeness, are shown in Figs.[T]and|2]respectively. 
In each figure, the top panel shows the mean flux expected, if 
the sources follow the large-scale structure distribution and the 
modulation pr discussed above is ignored. In the bottom panel, 
a modulation of (t)ioo - 1 has been adopted. Plotting this mean 
flux corresponds formally to what would be observed for a con- 
tinuous distribution of sources (i.e. infinite apparent source den- 
sity). The top panel, with no modulation, indicates what would 
be seen if the probability of scattering on the magnetized re- 
gions were everywhere unity, i.e. (t)i()o +00, or if the sources 
were steady emitters. The bottom panel shows the same map, 
but for bursting sources and including the modulation by the 
patchy distribution of magnetic field in the local Universe. For 
the sake of clarity, we have constructed these plots with a rela- 
tively narrow gaussian filter in HEALPix (Gorski et al. 2005) in 
order to smooth the galaxy distribution to a relatively high reso- 
lution ~ 1°. This resolution is higher than that of the PSCz cat- 
alog (about 7°) and smaller than the typical deflection expected 
for the interaction with magnetic fields. In the actual simulations 
discussed below, we used a degraded filter with 4° x 4° pixels. 

To take into account the effect of the magnetic configuration 
on the received flux, we multiply the mean flux from each indi- 
vidual cell (in three-dimensional Galactic coordinates r, I, b) by 
the probability pj. We use several values of (t)ioo: 0.1, 1, 3, and 
00. We then obtain different models of CR distribution depending 
on the optical depth value. In the case (t) 00, the probability 
Pj is equal to 1 and we are able to observe all p particles roduced 
from powerful bursting sources that follow the large-scale struc- 
ture (LSS). Hereafter, we refer to this configuration as "the LSS 
model". The two other magnetic configurations will be named 
according to the value of the optical depth, e.g. (t) = 3 model 
and (r) - 1 model. 

Figures [T] and 12] present the sky maps of the UHECR differ- 
ential flux normalized over the total expected flux for UHECR 
energy thresholds £ = 60 EeV and £ = 80 EeV. Generally 





Fig. 1. Sky maps in Galactic coordinates of the UHECR dif- 
ferential flux indicated in Eq. (O for an energy threshold E - 
60 EeV corresponding to a depth of 200 Mpc. Top panel: LSS 
model, in which sources follow the matter distribution according 
to the PSCz survey and in which (t)ioo +<^, corresponding 
to the assumption that the probability of scattering on a magne- 
tized system is unity; bottom panel: same, but with (r)i()o - 1. 
The grey zone represents the mask of the PSCz catalog in which 
there is no information. 



speaking, the distribution of the flux reflects the galaxy distri- 
bution following the LSS, but in the bottom map (correspond- 
ing to (r) = 1) the flux seems to be lower than in the top 
panel (corresponding to the LSS model, {t)ioo — » °°)', in the 
bottom panel, the smaller number of scattering centers in the 
voids concentrate the flux around regions of overdensity such 

as the Virgo cluster (b ~ 74°, I 80°), the Shapley cluster 

(b ~ 29°, I ~ -54°), the Centaurus cluster (b ~ 21°, I ~ -58°), 

the Hydra cluster (b ~ 26°, 1 101°), the Pavo-Indus cluster 

(b ~ -23°, / ~ -28°), the Fornax cluster {b ~ -53°, / ~ -124°), 
and the Perseus Pices cluster {b 17°, / ~ 124°), with a con- 
comittant reduction in the underdense areas of the sky. The gen- 
eral effect of the distortion that we discuss here is thus to cluster 
the arrival directions more strongly in the densest areas of the 
sky. 

3.2. Statistical tests 

To quantify the observed effect on the sky maps, we use vari- 
ous statistical tests that have been discussed in the literature and 
are related to the anisotropy of ultr ahigh energy cosmic ray ar- 
rival directions: the test (Cu ocol 20071). the one-d imensional 
Kolmogorov-Smirnov test (iKoers & Tinvakovll2009l) . the two- 
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For the Xc test, we divide the sky into angular bins of size 
4°x4° to avoid the local magnetic field effect as mentioned above 
and use two reference models: the isotropic model in which the 
arrival directions are purely isotropic (in the absence of experi- 
mental aperture sensitivity) and the LSS model with (t)ioo c»- 
We then compare the number of simulated events per bin in the 
tested model to average numbers of events from the two refer- 
ence models following the formula 



Xr 



(NJ - (M,LSs))((jV,',iso) - (M,LSS)) 
<M',LSS> 



(4) 




Fig. 2. Same as Fig. \T\ for an energy threshold E 
corresponding to a depth of 100 Mpc. 



80 EeV, 



dimensional Kolmogorov-Smirnov test dHarari et al ] I2009h . 
the Y test dKoers & Tinvakovl l2009h . the correlation test Xc 
jKashti & WaxmanI l2008h . and the Kuiper test. Although we 
tested them all, we present in what follows the results of those 
tests that have proven to be the most sensitive: the Xc test, the 
two-dimensional Kolmog orov-Smirnov test, and the Kuiper test. 
The Xc test was found in iKashti & WaxmanI (l2008h to be more 
sensitive to the anisotropy signature than the power spectrum 
and the two-point correlation function. The two-dimensional 
Kolmogorov-Smirnov test was used in Harari et al. (2009) to ex- 
clude isotropy when studying the UHECR distribution. And fi- 
nally the Kuiper test, wh ich is similar to the Y test defined in 
iKoers & Tinvakovl (l2009h . was found to have the highest statis- 
tical power of all considered statistical tests when considering 
a small number of high energy events. We tried to span dif- 
ferent types of tests: the Xc test relies on binning the sky, and 
the Kuiper and two-dimensional Kolmogorov-Smirnov tests de- 
pend on cumulative functions of the flux, the first being one 
dimensional, the second involving the Galactic coordinates I 
and b. Using these tests, we compare the distributions of simu- 
lated events from different models to the average distributions of 
events in the LSS model, (i.e. (r)ioo —> co)- 

As explained in Section 12.31 we generate two sets of sim- 
ulated events for each model, corresponding to two possible 
threshold energies E = 60 EeV and E - 80 EeV, i.e. a set of 
A^ev = 100 particles sampled neglecting the contribution from 
sources beyond 100 Mpc (corresponding to the horizon for par- 
ticles above 80 EeV) and a set of A^ev = 1000 particles from pos- 
sible sources within 200 Mpc (corresponding to the horizon for 
particles above 60 EeV). 



where NJ denotes the number of simulated events in angular bin 
; from the model to be tested with magnetic optical depth r, 
(M'.LSs) denotes the average number of events expected in an- 
gular bin / in the reference LSS model and (Mjso) corresponds 
similarly to the average number of events expected in angular 
bin / in the isotropic model. 

For the other two tests, we take as a reference model the LSS 
model and use cumulative distributions of the flux. The Kuiper 
test measures the quantity V, which is defined as the maximum 
deviation above an d below the two considere d cumulative distri- 
butions. Following iKoers & Tinvakovl (l2009l) . we name Crand(lF) 
the cumulative distribution of the simulated integrated flux from 
the model to be tested and Ctotai(!^) the cumulative distribution 
of the expected integrated flux from the reference model. We 
then have 

V = max[C,-and(r) - Ctotai(^n] + max[Ctotai(r) - Cra„d(n] -(5) 

For the two-dimensional Kolmogorov-Smirnov test, each 
simulated event n is characterized by its galactic coordinates 
{l„, b„y, the sky plane is then divided into four quadrants close 
to this direction (I < /„, b < bn), (I > l,,, b < b„), (I < l„,b > b,,), 
and (/ > l„,b > b„) and the test computes the differences be- 
tween the cumulative distributions of the tested and reference 
models in each of the four quadrants. This procedure is repeated 
4 X Nev times and D^, is defined as the maximum of these dif- 
ferences (between the cumulative distributions of the tested and 
reference models). In other words. 



Dks = max[|Ce (/,fe)-Cl,(/,fo)|]. 



e=i..4 



(6) 



where (l,b) are the Galactic coordinates of the cosmic ray 
positions and Q is one of the four quadrants with respect to a 
simulated event position in which data can be accumulated. 

With that being defined, we plot the histograms of the Xc, V, 
and Dks distributions from 1000 Monte Carlo realizations. In 
detail. Fig. |3] presents the distributions of Xc from 1000 Monte 
Carlo realisations of 100 simulated events of energy above 80 
EeV (top panel) and 1000 simulated events with energy above 
60 EeV (bottom panel). In the top panel, we note some over- 
lap between the LSS and the other models with finite (t)ioo, the 
magnitude of which depends on the value of <t)ioo- However, the 
effect of the finite magnetic optical depth becomes pronounced 
in the lower panel, representing the case of A^ev = 1000 above 
60 EeV. We note that these statistics are for instance well within 
the capabilities of the JEM-EUSO mission. In Fig. |4] we plot 
the distributions of V (top) and Dks (bottom) for 1000 Monte 
Carlo realisations of 1000 simulated events with energy above 60 
EeV; we have chosen this set of events for illustration because 
it appears more sensitive, as indicated by the previous figure. 
Qualitatively speaking, the Xc distribution appears to be a more 
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Fig. 3. Distribution of Xc from 1000 Monte Carlo realisations of 
100 simulated events with energy above 80 EeV (top) and 1000 
simulated events with energy above 60 EeV (bottom). Mean 
source density «s = lO^^Mpc"^. The dotted curves represent 
the isotropic model, the hatched curves denote the LSS model 
and from left to right we have {t)ioo = 0.1, (t)ioo = 1, and 
(r)ioo = 3. 



efficient discriminator of the effect discussed here. One might re- 
sort to definite statistical tests to quantify the sensitivity of these 
tests to the effect of finite optical depth but, given the exploratory 
nature of this discussion, the general trend indicated by the above 
figures suffices. 

The comparison of the top and bottom panels of Fig. [3]indi- 
cates that a larger number of events is more important in probing 
this effect, than a stronger anisotropy of the arrival directions. 

Figure|5]allows us to probe the sensitivity of the above results 
to the mean source density. In this figure, we show the distribu- 
tion of Xc from 1000 Monte Carlo realisations of 1000 simulated 
events of energy above 60 EeV for «s - 10"''Mpc"^, which cor- 



200 



150 



100 



50 




0.25 




0.45 



Fig. 4. Distribution of V (top) and Dks (bottom) from 1000 
Monte Carlo realisation 1000 simulated events with energy 
above 60 EeV. Mean source density «s = lO^^Mpc"^. The dot- 
ted curves represent the isotropic model, the hatched curves de- 
note the LSS model and from right to left we have {t)io() = 0.1, 
<t)i()() = 1, and <t)ioo = 3. 



responds to the existing lower bound on «s. When comparing 
with the bottom panel of figure [3] where - lO^^'Mpc"^, we 
note that a lower source density implies a relatively weaker sen- 
sitivity to the effect of finite magnetic optical depth. This results 
directly from the stronger clustering of events in given locations 
of the sky, which limits the number of directions with which 
one can sample the sky maps shown in Fig. [T] Several remarks 
are in order in this respect. First of all, one may hope that fu- 
ture large-scale experiments will be able to determine the mean 
source number density from the statistics of multiplets. If not, at 
the very least, one would infer that «s S 10"^ Mpc"^, hence for 
the present purpose the exact value of the present parameter is 
unimportant. Assuming that the mean source density is known, 
it would be interesting to explore the information contained in 
the statistics of multiplets in different regions of the sky, as these 
would obviously provide additional tools that are sensitive to the 
above effect; this a study requires us to handle a larger dimen- 
sionality of parameter space, thus is left to future work. 
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Fig. 5. Distribution of Xc from 1000 Monte Carlo realisations 
of of 1000 simulated events with energy above 60 EeV. The 
adopted mean source density is «s = lO^^Mpc"^. The dotted 
curves represent the isotropic model, the hatched curves denote 
the LSS model and from left to right we have (t)ioo = 0.1, 
<t)ioo = 1, and (t)ioo = 3. 



3.3. Application to tine Pierre Auger Observatory data 

To date, the only catalog of events above 60 EeV publicly 
available is that releas ed by the Pierre Auger Collaboration 
([Abraham et al.l l2008bh . We thus now use this catalog as a 
testbed for the previous tests. Since the catalog contains 27 
events above 57 EeV, which is reduced further to 18 by apply- 
ing the PSCz mask, the statistics is obviously limited and the 
present exercise should be taken as an example. 
We use the Auger aperture given by 



W{5) - cos(flo) cos(5) sin(a,„) + a,„ sin(flo) sin(5), 



where oq 
given by 



and 



-35° is the Auger southern site latitude and a,„ is 
[ Oif^>l 
I CO ' ) otherwise 



cosCflmax) - sin(ao) sin(g) 

COS(fl()) cos(5) 



Figure (|6|l shows the statistical distributions of Xc recon- 
structed for the different models using Monte Carlo simulations, 
while the vertical line indicates the value obtained for the Auger 
events. In accordance with the choice of threshold energy, we 
have limited the horizon to 200 Mpc. There is an almost com- 
plete overlap between the histograms of the LSS and the other 
models. Our study confirms that the statistics of these published 
Pierre Auger events appear inconsistent with the isotropic model 
and consistent with the LSS model (see Kashti & Waxman 200Ji 
iKoers & Tinvakovll2009l: lAublin et al.ll2009[). We also reco ver a 
(marginal) tendency, as noted in Kashti & WaxmanI (l2008l) . for 



the arrival directions to favor a scenario in which the appar- 
ent arrival directions cluster more strongly than the LSS traced 
through the PSCz; in our case, this is represented by the Xc test 
for the Auger events tending to measure smaller (t)ioo values, 
although we note that this tendency remains at most marginal. 
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Fig. 6. Distribution of Xc from 1000 Monte Carlo realizations 
of 18 simulated events of energy above 60 EeV. The mean 
source density is «s = lO^^Mpc"^. The dotted curves represent 
the isotropic model, the hatched curves denote the LSS model, 
and from left to right we have (r)i(K) - 0.1, (t)i()o - 1, and 
{t")ioo = 3. The vertical line corresponds to the value obtained 
for the Pierre Auger data. 

We noted in figures [3] and |4] that the studied effect becomes 
noticeable when the number of simulated events with energy 
above 60 EeV is of the order of 10^ - 10^; as expected, the effect 
is too weak to affect the Pierre Auger results and allow any clear 
distinction between the different models. 

4. Conclusions 

We have discussed the possible distortion of the sky maps of 
UHECR arrival directions associated with the inhomogeneous 
distribution of large-scale magnetic fields. As the effect of the 
extragalactic magnetic field on the UHECR trajectories may be 
described by a set of stochastic interactions with magnetized 
structures, its influence can be conveniently characterized by 
the optical depth to magnetic deflection r. Assuming that the 
UHECR sources are of the bursting type, such as gamma-ray 
bursts or newly born magnetars, and that their occurrence rate in 
the GZK sphere is low, we have argued that the relative lack of 
magnetized structures in those areas of the sky that have t < 1 
should lead to a depleted cosmic -ray flux from these directions 
compared to those detected from these directions if one assumed 
that the sources were steady emitters. This is mainly related to 
the particles in those areas of the sky, having a probability t of 
interacting with a magnetized structure, which in the above sce- 
nario represents a necessary condition for the source to be ob- 
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servable in an experiment lifetime. It is important to emphasize 
that the effect that we have discussed here is different from the 
standard angular deflection associated with extragalactic mag- 
netic fields. While angular deflection modifies the arrival direc- 
tion without altering the flux, the effect that we discuss directly 
affects the flux from some regions of the sky. 

To quantify this distortion, we have constructed sky maps of 
UHECR arrival directions from bursting sources residing in the 
large scale structure for two energy thresholds of 60 EeV and 
80 EeV (see Figs. [1] and |2]i. We have considered different config- 
urations of the large-scale magnetic field and taken into account 
the modulation associated with the probability of experiencing at 
least one interaction with a magnetized system = 1 -exp(-T). 
Each model is then characterized by the value of (t)ioo, which 
expresses the average optical depth of magnetic deflection cal- 
culated to a depth of 100 Mpc. We have then used various statis- 
tical tests previously described in the literature, particularly the 
Xc test [see Eq. ©I, which measures the correlation between 
the simulated events from the model to be tested and the aver- 
age expected number of events from two reference models, one 
model with isotropic arrival directions and one model in which 
the sources are distributed according to the large-scale structure 
such that the effect of finite magnetic optical depth is neglected. 
Inspired by the statistics expected for future large-scale exper- 
iments, we have applied these tests to two sets of events: A^ev 
= 100 events with £ > 80 EeV and A^ev = 1000 events with 
£ > 60 EeV, respectively. The histograms from the second set 
of events show more distinguishable models because of the high 
number of expected events (see Figs.[3]and|4]i. We have also con- 
sidered the catalog of events (27 events above 57 EeV) released 
by the Pierre Auger Observatory as an example. 

In conclusion, we have found that the distortion examined 
here cannot affect present scale experiments but should be con- 
sidered when performing anisotropy studies with future large- 
scale experiments. We have noted a slight sensitivity of the dis- 
tortion strength to the mean ultrahigh energy cosmic ray source 
density, which would however disappear if this latter were deter- 
mined using the statistics of multiplets. 
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